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Abstract. We have performed numerical simulations of inertial particles in random 
model flows in the white-noise limit (at zero Kubo number, Ku = 0) and at finite Kubo 
numbers. Our results for the moments of relative inertial particle velocities are in 
good agreement with recent theoretical results (Gustavsson & Mehlig 2011a) based on 
the formation of phase-space singularities in the inertial particle dynamics (caustics). 
We discuss the relation between three recent approaches describing the dynamics and 
spatial distribution of inertial particles suspended in turbulent flows: caustic formation, 
real-space singularities of the deformation tensor, and random uncorrelated motion. 
We discuss how the phase- and real-space singularities are related. Their formation 
is well understood in terms of a local theory. We discuss implications for random 
uncorrelated motion. 
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1. Introduction 

The dynamics of particles suspended in randomly mixing or turbulent flows ('turbulent 
aerosols') has been studied intensively for several decades. Recently, substantial progress 
in understanding the dynamics of turbulent aerosols has been achieved (see the papers 
published in this special issue and the references cited therein). 

The phenomenon of spatial clustering of independent point particles subject 
to Stokes drag in turbulent flows is now well understood: below the dissipative 
length scale (where the fluid flow is smooth) the particles eventually cluster onto 
a fractal set in configuration space. The corresponding fractal dimension has been 
determined by means of direct numerical simulations (Bee 2003) as well as theoretical 
approaches (Wilkinson et al. 2007, Gustavsson & Mehlig 20116). Different mechanisms 
['preferential concentration' (Maxey 1987) and 'multiplicative amplification' (Wilkinson 
et al. 2007, Gustavsson & Mehlig 20116)] contribute to spatial clustering. A third 
mechanism giving rise to particle clustering was recently studied by following the 
deformation of an infinitesimally small volume of particles transported along a particle 
trajectory ['full Lagrangian method' (Ijzermans et al. 2010)]. The small volume may 
vanish at isolated singular points in time, giving rise to instantaneous singularities in 
the particle-concentration field. Using this approach the statistical properties of these 
singularities were analysed by Meneguz & Reeks (2011). 

One important reason for studying spatial clustering of inertial particles is that this 
phenomenon is argued to enhance the rate at which collisions occur in turbulent aerosols 
at small values of the 'Stokes number'. This dimensionless parameter, St = (77-) _1 , is 
given in terms of the particle damping rate 7 and the relevant correlation time r of the 
flow. Both are defined below. 

Arguably spatial clustering has an effect upon the collision rate at small Stokes 
numbers. But there is a second mechanism that leads to a significant enhancement 
of the collision rate as the Stokes number increases: direct numerical simulations of 
particles suspended in turbulent flows (Sundaram & Collins 1997, Wang et al. 2000) show 
that relative particle velocities at small separations increase substantially as the Stokes 
number is varied beyond a threshold of order unity. In (Falkovich et al. 2002, Wilkinson 
et al. 2006) this behaviour was explained by the occurence of singularities in the particle 
dynamics, causing large relative velocities at small separations. These singularities 
occur as the phase-space manifold folds, as illustrated in Fig. [TJ As a consequence, 
particle velocities at a given point in space become multi-valued, causing large velocity 
differences between nearby particles. The boundaries of the folding region are referred 
to as 'caustics' (Wilkinson et al. 2005, Crisanti et al. 1992). It was shown that 
the rate of caustic formation is an activated process (Wilkinson et al. 2005, Duncan 
et al. 2005, Gustavsson & Mehlig 2012). This explains the sensitive dependence of 
the rate of caustic formation upon the Stokes number observed in direct numerical 
simulations of particles in turbulence (Pumir & Falkovich 2007). 

An alternative way of characterising relative velocities of inertial particles was 
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Figure 1. a Trajectories of a one-dimensional model, particle positions as a function 
of time. Also shown (b, c, d): phase-space manifolds (velocity v versus position x) 
demonstrating how the phase-space manifold folds over at a caustic. Panels a to d are 
similar to Fig. 1 in Gustavsson & Mehlig (2011a). Also shown (e, f, g): position x as 
a function of initial position xq. Parameters: St = 300, Ku = 0.1. 

suggested in (Fevrier et al. 2005, Simonin et al. 2006). The authors of these papers 
decomposed inertial particle velocities into two contributions: a spatially correlated, 
smoothly varying 'filtered' velocity field, and a random, spatially and temporally 
uncorrelated contribution, commonly referred to as 'random uncorrelated motion' 
(Reeks et al. 2006, Masi et al. 2011). 

The aim of this paper is twofold. First we summarise results of numerical 
simulations of particles suspended in model flows (Figs. |3]- 18]). Our numerical results 
for the moments of relative velocities of inertial particles are in quantitative agreement 
with recent analytical results based on the notion of caustic formation (Gustavsson & 
Mehlig 2011a). Second we demonstrate that caustic formation not only provides an 
understanding of relative velocities at small separations, caustic formation also explains 
spatial clustering due to singularities in the local deformation tensor, and the existence 
and properties of random uncorrelated motion. 

We conclude the introduction by summarising our results in more detail. In this 
paper we show that recent predictions by Wilkinson et al. (2006) and Gustavson & 
Mehlig (2011a) based on the notion of caustic formation describe many aspects of the 
fluctuations of relative velocities at small separations. We compare formulae for the 
moments of relative velocities (Eqs. (fT8l) and (1201) below) to new results of numerical 
simulations of one- and two-dimensional models for inertial particles suspended in white- 
noise flows, and for a three-dimensional kinematic simulation of particles suspended in 
an incompressible flow field with an energy spectrum typical of the small scales of 
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Figure 2. (Left) multi-valued velocities of particles suspended in a two-dimensional 
random flow with finite Ku and St as described in Section EOl The base of each red 
arrow corresponds to a particle position (taken to be on a regular grid in the x-y- 
plane). The orientation of the velocity is that of the arrow. All arrows have the same 
length, the magnitudes of the velocities are not shown. The blue line delineates the 
position of the caustics in the x-j/-plane. The region of multi-valued velocities ends 
in a cusp that is only approximately resolved. In Section |4] it is explained how multi- 
valued velocities between caustics give rise to so-called random uncorrelated motion. 
Parameters: Ku = 1, St = 10. (Right) particle-density in the x-y- plane, showing 
significantly enhanced particle-number density in the vicinity of the caustic line. Same 
parameters as above. Black corresponds to high density, white to low density. 



turbulence. We find good agreement. This demonstrates that Eqs. ([181) and (|20|) which 
were derived in the white- noise limit, are valid more generally. 

Further we examine the prediction by (Fevrier et al. 2005, Simonin et al. 2006) that 
the so-called longitudinal second-order structure function of relative velocities tends to 
a finite value at vanishing separations in the presence of random uncorrelated motion. 
The analytical theory, Eqs. ffl~8]) and ( 120]) below, shows that this is true for sufficiently 
large Stokes numbers [the case examined numerically by (Simonin et al. 2006)]. But 
at Stokes numbers smaller than a critical value, the structure function tends to zero, 
despite the fact that there may still be a substantial singular (multi-valued) contribution 
to relative velocities due to the formation of caustics. 

We discuss in detail that the singularities of the deformation tensor are in fact 
caustic singularities, as pointed out by (Wilkinson et al. 2007). We study the dynamics 
of the deformation tensor JJ, and the matrix Z of particle-velocity gradients. We show 
that as det J approaches zero, TrZ — > — oo. We briefly remark upon the statistical 
properties of the singularities (Meneguz & Reeks 2011). 

In summary, we demonstrate that the notion of random uncorrelated motion, and 
the occurrence of zeroes in the local deformation tensor can both be explained in terms 
of caustic formation, both qualitatively and in many ways quantitatively. Last but not 
least our results indicate that the white-noise approximation successfully describes many 
aspects of turbulent aerosols. 
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The remainder of this paper is organised as follows. In Section [2] we introduce 
the models analysed in this paper: inertial particles suspended in a two-dimensional 
incompressible random flow in the white- noise limit, and a kinematic simulation of 
inertial particle dynamics. Section [3] summarises what is known about the rate of caustic 
formation and discusses consequences for the fluctuations of relative particle velocities. 
We compare the analytical theory to results of numerical simulations of the models 
described in Section |2j In Section H] we briefly review the notion of random uncorrelated 
motion, and compare the conclusions of (Fevrier et al. 2005, Simonin et al. 2006) to 
our analytical and numerical results. In Section [5] we describe the dynamics of the local 
deformation tensor and its correspondence to the dynamics of the matrix of particle- 
velocity gradients. Finally, Section [6] contains our conclusions. 

2. Model 

The motion of small, non-interacting spherical particles suspended in a flow is commonly 
approximated by 

r = v , v = 7(14 — v) . (1) 

Here r and v are the position and velocity of a particle, u(r,t) is the velocity field 
evaluated at the particle position, 7 is the viscous damping rate, and dots denote time 
derivatives. The components of the vector r are denoted by rj, j = l,...,d in d 
dimensions. The components of u and v are referred to in an analogous way. Sometimes 
it is more convenient to denote the components of r by (x,y,z) instead of {ri,r 2 ,r 3 ). 
We use the two notations interchangeably. 

For Eq. (CQ) to be valid, it is assumed that the particle Reynolds number is small, 
that Brownian diffusion of the particles can be neglected, and that the particle density 
is much larger than that of the fluid. We also assume that the velocity field u varies 
smoothly on small spatial and temporal scales with smallest length-and time scales rj 
and r (the Kolmogorov scales for turbulent flows). The typical magnitude of the velocity 
field is denoted by u . 

In dimensionless units (t = t'/^f, r = rjr', v = 7771/, u = 'jrju' and dropping the 
primes), the equation of motion becomes 

r = v , v = u — v . (2) 

The Stokes number does not appear explicitly in this equation, but the fluctuations 
of the dimensionless velocity u depend upon St (see Eq. (TjJ below). In addition to 
the Stokes number, the dynamics is characterised by a second dimensionless number, 
the 'Kubo number' Ku = uqt/7]. We note that turbulent flows have Ku ~ 1. In 
the remainder of this paper we frequently refer to these two dimensionless numbers. 
For a discussion of further dimensionless parameters see (Wilkinson et al. 2007). The 
numerical results shown in the following were obtained for two different models. These 
models are introduced in the following two subsections. 
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2.1. Random-flow model 

Following (Wilkinson & Mehlig 2003, Wilkinson et al. 2005, Duncan et al. 2005, 
Wilkinson et al. 2007) we approximate the incompressible velocity field u(r,t) in 
Eq. d2J) by a Gaussian random function that varies smoothly in space and time. We 
discuss results for one- and two-dimensional versions of the random-flow model. The 
one- dimensional case is most easily analysed, the two-dimensional incompressible case 
is important (since one-dimensional flows are special, they are always compressible 
which give rise to a path-coalescence transition (Wilkinson & Mehlig 2003)). A two- 
dimensional incompressible velocity field can be written in terms of a stream function 
if>(r,t): 

u(r,t) = V Atjj(r,t)e 3 . (3) 

Here e3 is the unit vector _L to the x-y-plane. We assume that ip(r,t) is a Gaussian 
random function with (ip) = and correlation function 

(^(r,t)V>(0,0)) = ^Ku 2 St 2 exp [-\r\ 2 /2 - St |t|] . (4) 

in dimensionless variables. 

In this paper we also refer to results of a one-dimensional random-flow model. This 
is defined in an analogous fashion in terms of a Gaussian random flow velocity u(x,t) 
with zero mean and correlation function 

(u(x, t)u(0, 0)) = Ku 2 St 2 exp [-x 2 /2 - St \t\] . (5) 

We note the one-dimensional flow is compressible. The numerical data shown in Figs. CD 
and |2] are obtained by computer simulations of the models described above. 

We simplify the model by linearising Eq. (j2J). This yields the following equation for 
the dynamics of a small separation R = r\ — r 2 and velocity difference V = V\ — v 2 
between two particles: 

R = V , V = AR - V . (6) 

Here A is the matrix of fluid velocity gradients, with elements = dui/drj. 

To simplify further, we take the white-noise limit of this model. This limit 
corresponds to 

Ku — > and St — > oo such that e 2 = QKu 2 St = const. (7) 

Here e is a dimensionless measure of the particle inertia introduced by (Mehlig & 
Wilkinson 2004) [see also (Wilkinson et al. 2007)]. We take c\ — \ for one-dimensional 
flows [this is consistent with the convention used in (Gustavsson & Mehlig 2011a)]. 
For incompressible two-dimensional flows we take c 2 = 1/2, as in (Gustavsson & 
Mehlig 20116). In the white- noise limit, the instantaneous value of the velocity gradient 
A in (|6]) becomes independent of the particle position. In two spatial dimensions, we 
denote the independent random increments of the elements An, A\ 2 and A 2 \ of A in 
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The results shown in Figs. El HI and |5] are obtained by computer simulations of this 
model, approximating the time-dependence of A(r(t),t) as a white-noise signal. 

2.2. Kinematic simulation 

As an alternative to the single-scale white-noise model introduced in the previous 
subsection, we simulate a turbulent incompressible velocity field in a three-dimensional 
periodic box by a large number of Fourier modes varying randomly in space and time. 
The modes are chosen in such a way that the associated energy spectrum approximates 
a prescribed form, namely that originally used by (Kraichnan 1970). The model is 
identical to that used by (Ijzermans et al. 2010, Meneguz & Reeks 2011). For convenience 
we briefly summarise its relevant features below. For details, we refer the reader to 
(Ijzermans et al. 2010, Meneguz & Reeks 2011). 

In dimensionless form, the incompressible velocity field u(r,t) is represented as a 
Fourier series of ./V modes (N = 200 in our simulations): 
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with random coefficients and b^ n \ random wave numbers k^ n \ and random 

frequencies oj^ n \ In order to guarantee the periodicity of the flow in a cube of dimensions 
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As mentioned above, the energy spectrum E(k) is taken to be (Kraichnan 1970): 
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This spectrum is representative for low-Reynolds-number turbulence (Spelt & 
Biesheuvel 1997). The maximum of E(k) is located at k — 1 and the total kinetic 
energy J °° E{k)dk = 3/2. This corresponds to 3mq/2 in dimensional form. The use 
of the Kraichnan energy spectrum results in a relatively small separation of scales; in 
our simulations, the smallest wave number fcW ~ 0.25 and the largest wave number 
^ 2.14. The frequencies u( n ' are chosen randomly from a Gaussian distribution 
with zero mean and a variance proportional to k^ 1 '. This implies that the Kubo number 
is of order unity. Following (Spelt & Biesheuvel 1997), we take the variance to be0.4fcM. 
Finally, the coefficients and are determined by choosing a random direction in 
Cartesian space, and by picking a length randomly from a Gaussian distribution with 
zero mean and a variance 9/(2N). By doing so, the mean kinetic energy at a given 
position in space 

E km {r) = ]-Yim \ [ dt\u(r,t)\ 2 (14) 

N 1 r i 
= V — ^- loW Afc (n) | 2 + |b W Afc (n) | 2 , 

n=l I I 

is approximately equal to 3/2 for all values of r. 
3. Caustics 

3.1. One spatial dimension 

As illustrated in Fig. HJ caustics form when the phase-space manifold folds over. In one 
spatial dimension this happens when the slope of the manifold becomes infinite, that is 
when z = dv/dx — > —oo. The rate at which this occurs is determined by the equation 
of motion for z (Wilkinson & Mehlig 2003): 

z = A-z-z 2 . (15) 

Here A = du/dx represents the random driving by the fluid- velocity gradients. In the 
case of independent particles (which we consider here), z goes through infinity in a 
symmetrical fashion. At large values of |z|, the random driving can be neglected, so 
that z ~ — z — z 2 . The corresponding deterministic probability distribution of z reads 
p(z) = C/[z(l + z)], and is valid in the tails of z. 

In the white- noise limit, Eq. ( Fl5l) is equivalent to a Fokker-Planck equation for the 
distribution of z. In (Wilkinson & Mehlig 2003) this equation was solved in one spatial 
dimension. The resulting rate of caustic formation [called 'rate of crossing caustics' by 
Wilkinson & Mehlig (2003)] can be written as (Gustavsson & Mehlig 2012): 



J caustic 1 f-j- 

= — lm 

7 2n 
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where e 2 = Ku 2 St (see Section |2~T|) . In Eq. (1161) . Ai(y) is the Airy function. In the limit 
of small values of e, this expression exhibits the asymptotic behaviour 

Jcaustic _ 1 _ c~ 1/(6e2) . (17) 



7 V2tt 

Eq. (fTB"l) shows that the number of caustics increases rapidly as e 2 passes through 1/6 
(Wilkinson et al. 2005, Wilkinson et al. 2006). This sensitive dependence is commonly 
referred to as an 'activated law', in analogy with the sensitive temperature depdence of 
chemcial reaction rates in Arrhenius' law. Gustavsson & Mehlig (2012) computed the 
one-dimensional rate of caustic formation at small but finite Kubo numbers and found it 
to sensitively depend on the Stokes number: in this case too, the St-dependence exhibits 
'activated form': J ca ustic/7 ~ exp[— S'(St)/Ku 2 ], where S is an St-dependent 'action'. In 
the white-noise limit, S = l/(6St), consistent with Eq. ffTTl) . 

As Fig. [1] shows, particle- velocities become multi- valued between two caustics in the 
wake of a singularity, giving rise to large relative velocities between nearby particles. 
While the rate of caustic formation is determined by the rate at which the local 
quantity z = dv/dx tends to — oo, the distribution of relative velocities at small particle 
separations is determined by the solution of the full non-local equations ([6]) for particle 
separations and relative velocities (Gustavsson & Mehlig 2011a). 

A consequence of large relative velocities at small separations is that between 
caustics, particles collide frequently with large relative velocities (c.f. Fig. [1]), giving 
rise to a large collision rate (we note, however, that in this paper it is assumed that the 
particles are independent point particles that do not actually collide). 

By contrast, in the absence of caustics, particles may still approach each other due 
to fluctuations of the underlying flow-velocity field. At small separations the flow is 
smooth, and in this regime relative velocities between particles are expected to tend to 
zero as the particles in question approach each other. 

Which one of these two mechanisms of bringing particles together makes the 
dominant contribution to the collision rate depends upon the value of St and on the 
particle size a (separation 2a at the point of contact). Relative velocities of particles 
thrown at each other due to the formation of caustics are expected to make the 
dominant contribution to the collision rate if St is large and/or when the particles 
are sufficiently small. Particles slowly approaching each other ('logarithmic diffusion') 
dominate otherwise. 

In the white- noise limit, and in one spatial dimension, an asymptotic approximation 
for the moments of relative velocities at small separations was derived in (Gustavsson 
& Mehlig 2011a): 

/oo 
dV\V\ p p(X, V) ~ i^Xr^- 1 + C p . (18) 
-oo 

Here X = x\ — x 2 and V = V\ — v 2 are the separation and the relative velocity of 
a pair of particles, and p(X, V) is their distribution function. It is assumed that 
IX I <§C 1 and p > —1. Further, D 2 is the correlation dimension of the phase-space 
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Figure 3. Moments of the radial velocity m p (R) plotted against distance R for two 
different values of e: e = 0.03 (left) and e = 0.06 (right). Data from numerical 
simulations of the two-dimensional white-noise model described in Section 12.11 are 
shown as markers. The correlation dimension c?2 and the coefficients B p and C p in 
the small R approximation (|20[) are numerically fitted to the data in the interval 
bounded by vertical black dashed lines. The resulting moments for small R ([20)) are 
shown as solid lines. The caustic contribution C p i? d_1 (dashed dotted) and the smooth 
contribution B p R p+d2 ^ 1 (dashed) are also shown. Parameters: p = (red o), p = 1 
(green □), p = 2 (blue 0) and p = 3 (magenta A). 

attractor and B p and C p are mo del- dependent constants [which were not derived by 
(Gustavsson & Mehlig 2011a)]. The form of Eq. ([IB]) is consistent with the form inferred 
from simulations of relative-particle dynamics in a one-dimensional Kraichnan model 
(Cencini 2009). 

The second term in Eq. f [T8|) . C p , is due to multi- valued velocities between caustics. 
This contribution, in one spatial dimension, does not depend upon \X\ for small values 
of \X\. In other words: it remains finite as \X\ — > 0. This is a consequence of the fact 
that as the manifold in Fig. [TJ folds over, particles initially far apart are thrown at each 
other quickly. 

The first term in Eq. ( !T8l) . 5 P |X| P+D2_1 , vanishes as \X\ — >■ 0. It constitutes 
the main contribution to m p (X) in the absence of caustics and is affected by spatial 
clustering: for a given value of p, the exponent is smallest (and thus the contribution 
largest) when D-i attains its minimum as a function of Stokes number. 

In Eq. f[T8|) the case p — 1 is of particular importance, since mi(X) is closely related 
(yet not identical) to the collision rate between particles at small separations X = 2a. It 
is expected that the coefficient of the caustic contribution in Eq. ([T8|) . C\, is proportional 
to the caustic formation rate J caU stic (Wilkinson et al. 2006). 

3.2. Two and three spatial dimensions 

In two and three spatial dimensions the caustic rate can be found in a way similar to the 
one-dimensional case (Wilkinson et al. 2007, Gustavsson & Mehlig 20116): the matrix 
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Z with elements Z^- = dvi/drj obeys the equation: 
Z = A - Z - Z 2 . 



(19) 



Here A is the matrix of fluid- velocity gradients introduced in Section El with elements 
Aij = dui/drj. In analogy with the one-dimensional case, tr(Z) — > — oo as caustics 
are formed. In the white- noise limit, we expect that the rate of caustic formation 
is again given by i fTTj) . In (Wilkinson et al. 2005, Duncan et al. 2005, Wilkinson 
et al. 2007) numerical factors in Eq. (|T7j) slightly different from 1/6 were quoted in 
two and three spatial dimensions. More recent numerical results (not shown) show that 
the asymptote (TTTI) is approached very slowly as e becomes small. Our best estimates 
at the smallest values of e indicate that the factor in the argument of the exponential in 
(fT71) is asymptotically the same (equal to 1/6) in one, two, and three spatial dimensions. 

Moments of relative velocities in two and three spatial dimensions obey laws 
analogous to ( |T8l) . At small separations (R <C 1) Gustavsson & Mehlig (2011a) found 



Here R — \R\ and vr = V-e^ is the radial projection of the relative velocity between two 
particles at separation R. Further, d 2 is the spatial correlation dimension, it is assumed 
that the Stokes number is small enough so that d 2 < d. As in the one- dimensional result, 
Eq. (Tl8|) . there are two contributions to the moments of relative velocities [compare the 
parameterisation of the St-dependence of the collision rate suggested by Wilkinson et 



The second term in Eq. ( 1201) is due to multi-valued velocities between caustics. But 
note that in two and three spatial dimensions, not all particle pairs thrown together give 
rise to close approaches. The reason is that in addition to having one relative coordinate 
pass zero at finite relative velocity (so that a caustic occurs), the other coordinates must 
be small, i.e. only particles heading sufficiently towards each other as the caustic occurs 
end up at small enough separations to contribute to the small R velocity moments. 
This explains the geometrical factor R d ~ l in (120]) . It is absent in one spatial dimension, 



Fig. |3] shows comparisons of Eq. ( 120]) with results of numerical simulations of the 
random-flow model described in Section I2TT1 The parameter d 2 in Eq. (120]) is determined 
as follows. Setting p = in ( 120]) and taking the limit R — > defines the spatial 
correlation dimension d 2 . The latter is found numerically by fitting mo(-R) to the power 
law R^- 1 . 

We now describe how the fits in Fig. [3] were obtained. The parameter d 2 was taken 
from Fig. HI The coefficients B p and C p in Eq. (120]) were fitted to the numerical results 
for different parameter values. The fitting region (the range of R over which Eq. ( 120] ) 
is fitted) lies between the dashed lines in Fig. [3j We observe good agreement between 
the numerical results and fits to Eq. (1201) . In particular, the results clearly show that 
the moments m p scale as R d ~ l for small values of R, independently of p. Fig. [5] shows 
the coefficient C\ of the caustic contribution obtained in this way as a function of e~ 2 . 




(20) 



al. (2006)]. 



d = l. 
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Figure 4. Spatial correlation dimension d,2 (°) as a function of e 2 for the model 
described in Section 12.11 The dashed red line shows the small-e theory discussed in 
(Wilkinson et al. 2010, Bee et al. 2008). 
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Figure 5. Amplitude C\ from numerical fits to m,\ in (|2"0")) (symbols) as a function 
of er 2 . Numerical simulations of the model described in Section |2~T1 The asymptotic 
St-dependence of the rate of caustic formation, (Tl7|) . is shown as a dashed line. 



Since this contribution requires the formation of caustics, we expect C\ to exhibit an 
e-dependence of the form (fT7]) . Fig. [5] shows that this is indeed the case. 

Fig. |6] shows results for m Q (R) and mi(R) obtained by kinematic simulations of the 
random-flow model described in Section l2~2"l for two values of the Stokes number, St = 0.4 
and St = 0.7. As expected, m (R) is of power-law form, reflecting spatial clustering. 
The corresponding correlation dimensions are shown, as a function of St, in Fig. [7J The 
correlation dimension exhibits the expected minimum (here at St ~ 0.4). Corresponding 
results for direct numerical simulations of particles in turbulent flows have been obtained 
by a number of authors (Wilkinson et al. 2010, Bee et al. 2010, Chun et al. 2005). 

The green squares in Fig. [6] correspond to numerical results for m>i(R) as a function 
of R. Consider first the left panel (St = 0.4). At small separations R we expect that 
mi(R) should scale as R d ~ l = R 2 , while it should scale as R d2 ~ R 2A at large values 
of R. Despite the fact that the two powers are rather similar, the two scalings can be 
distinguished in Fig. El In the right panel (St = 0.7), the caustic contribution R^ 1 
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Figure 6. Same as Fig- El for the model described in Section I2T21 with St = 0.4 (left) 
and St = 0.7 (right), and p = 0, 1. 
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Figure 7. Numerical results for the spatial correlation dimension 0*2 as a function of 
the Stokes number for the model described in Section l2~2l The dashed line shows a fit 
of the form d 2 = 3 - 12 St 2 . 



dominates. 

Given the data available from the kinematic simulations, it is more difficult to 
reliably determine the St-dependence of C\ by fitting (solid green line in Fig. |6]). Our 
best estimates are shown in Fig. |SJ The fits and the corresponding error bars were 
obtained by a non-linear least-squared fit using MATLAB 2011. We find that C\ 
depends very sensitively on St, as expected because the formation of caustics is an 
activated process. We expect (Gustavsson & Mehlig 2011a, Gustavsson & Mehlig 2012) 
that the St-dependence of C\ follows the law J caust i C /7 ~ exp[— S'(St)/Ku 2 ]. However, 
the range of Stokes numbers for which C\ can reliably be estimated is too small to 
determine the form of the function S'(St). 

Fig. El demonstrates that the magnitude of relative velocities at small separations 
depends very sensitively on the Stokes number. We argue that this is a consequence of 
the sensitive St-dependence of the rate of caustic formation. This explains the sensitive 
dependence on the Stokes number of collision velocities and collision rates of particles 
suspended in turbulent flows (Sundaram & Collins 1997, Wang et al. 2000, Zaichik & 
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Figure 8. Amplitude C\ from numerical fits to m\ in ([20]) as a function of 1/St 
according to numerical simulations of the model described in Subsec. 12.21 



Alipchenkov 2003). 

4. Random uncorrelated motion 

Singularities in the inertial-particle dynamics (corresponding to the formation of 
caustics) give rise to multi-valued particle velocities at locations in space bounded by 
caustics: any identical particles that are very close may move at substantially different 
velocities. Figs. [U and |2] illustrate this fact in one and two spatial dimensions. This 
implies in particular that the relative motion of inertial particles cannot be captured in 
terms of a 'hydro dynamic' approximation describing the particle velocities in terms of a 
smooth velocity field. In particular Fevrier et al. (2005) infer from their DNS calculations 
of inertial-particle motion in a homogeneous isotropic and stationary turbulent flow field 
that the velocity of a particle at a position r(t) at time t in a single realisation of the 
carrier flow field u(r, t) is given by the sum of two components 

v(r(t)=r,t) = v(r,t) + 5v(r(t) = r,t). (21) 

Here v(r, t) is a smoothly varying filtered velocity field which Fevrier et al. (2005) refer 
to as the 'mesoscopic Eulerian particle velocity field'. Values for the smooth component 
in any realisation are found by dividing the spatial domain into cells and calculating 
the average velocity associated with the number of particles in each individual cell 
(the number of particles in each cell being sufficiently large to form a statistically 
stationary average). The residual component Sv is termed the 'quasi Brownian velocity 
distribution component' by Fevrier et al. (2005). It is now commonly referred to as 
'random uncorrelated motion', or RUM for short (Reeks et al. 2006, Masi et al. 2011). 
This residual RUM part is assumed to be uncorrelated with the smooth part and with 
itself at infinitesimally small separations in space and time. 

The existence of multi-valued velocities between caustics is consistent with a 
singular contribution to the particle velocities, of the form of Eq. ( l2~TT) . We infer that 
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the extent of random uncorrelated motion [its relative contribution compared to the 
smooth part in Eq. (T2T]) ] must depend sensitively on the value of St, since the rate of 
caustic formation exhibits this sensitive dependence on the Stokes number. 

Let us consider the implications of Eq. (121]) and the accompanying assumptions 
for the second moment of the relative radial velocity between two particles, v R = 
(t>i - v 2 ) ■ e R : 

(«£) = + <*«i - «2 - Sv 2 ) ■ e R } 2 } 

= ([(v, - v x ) ■ e R ] 2 ) + ([5 Vl ■ e R ] 2 ) + ([5v 2 ■ e R ] 2 ) . (22) 

This result is of the same form as Eq. ( !20l) . The two right-most terms in (1221) correspond 
to the caustic contribution in (120]) . In other words, Eq. (1201) provides a quantitative 
prediction for the contribution of random uncorrelated motion to the moments of relative 
radial velocities. Consider for example the form of the so-called 'longitudinal structure 
functions' for relative velocities of the suspended particles. Simonin et al. (2006) argue 
that the second-order structure function remains finite as the spatial separation R 
between particle velocities tends to zero. In the notation of the previous section, the 
second-order structure function is given by 

^=3' < 23 > 

The limiting behaviour of s™(R) can be deduced from Eq. ( |20|) . From this equation 
we see that m (R) ~ ^nnn{d2,d}— 1_ The correlation dimension d 2 saturates to d at a 
critical Stokes number, St c (c.f. Fig. H] where d 2 = d for e 2 > e 2 ~ 1). For St > St c the 
suspended particles are uniformly distributed in space [see also (Bee et al. 2010, Salazar 
& Collins 2012)]. Let us consider this case. As R — > 0, the caustic contribution C 2 R d ~ 1 
to m 2 (R) dominates in Eq. (120]) . This implies that 

s {2 \R) -> const, as R , (24) 

as argued in (Simonin et al. 2006). For St < St c , by contrast, we find 

s {2) (R)^0 zsR^O. (25) 

More precisely, s ( - 2 '(R) tends to zero as g~ l (R) when R — > (the pair correlation function 
g{R) is given by g(R) = m {R)/R d - 1 ). 

We emphasise that the behaviour ( )24|) of the structure function in two and three 
spatial dimensions must be distinguished from the fact that the moments m p (X) of 
relative velocities in one spatial dimension always approach a positive constant as 
\X\ — > when St > 0. Indeed, we have shown that s^(R) may approach zero as 
R — > 0, yet multi- valued particle velocities do still give rise to a substantial singular 
contribution to the moments of relative velocities, as a consequence of singularities 
giving rise to caustics. 

Let us compare these findings to the results shown in Fig. 3(a) in Simonin et al. 
(2006). The data shown in this figure (except perhaps the data set labeled '1') imply 
that the structure function approaches a positive constant as R — > 0. We conclude that 
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the data sets shown (possibly with the exception of '1') correspond to Stokes numbers 
larger than St c . It should be noted that Simonin et al. (2006) define their Stokes number 
Sti in terms of the integral time scale of the turbulent flow. Here and in a large part of 
the literature on inertial particles in turbulent flows the Stokes number St is defined in 
terms of the Kolmogorov time r. Since usually St 3> Stz, it is plausible that most data 
sets in Fig. 3(a) in Simonin et al. (2006) correspond to St > St c . 

We conclude by noting that it has been shown (see (Mehlig et al. 2005, Wilkinson 
et al. 2007) and references cited therein) that the maximal Lyapunov exponent 
describing the dynamics of inertial particles suspended in incompressible flows is 
positive. This implies that the inertial particle dynamics is chaotic. In the limit of very 
large Stokes numbers, inertial particle dynamics is thus similar to the random motion 
of molecules in a gas [gas-kinetic limit, see (Abrahamson 1975)]. This justifies the view 
that there is a random uncorrelated component to the inertial particle dynamics. It is 
a consequence of the formation of caustics. 

5. Singularities in particle concentration 

Changes to the local concentration of inertial particles suspended in mixing flows can 
be described by the deformation tensor J with elements Jjj = dri/drj(0) evaluated 
along a particle trajectory r(t) with initial position r(0). The matrix J describes the 
relative motion of infinitesimally close particles. In particular, the volume spanned by 
the separation vectors between d+1 infinitesimally close particles in d spatial dimensions 
is given by 5V = \J\SVo, where J = det(J) and 5Vo is the initial volume, see Fig. [91 

Nothing prevents J from occasionally changing sign. This implies that the volume 
SV may shrink to zero, giving rise to a singularity in the local particle concentration 
oc 5V~ l (Wilkinson et al. 2005, Wilkinson et al. 2007, Ijzermans et al. 2010). The 
singularities influence the tails of the distribution of local particle concentration, making 
particle clustering highly non-Gaussian and intermittent (Meneguz & Reeks 2011). The 
zeroes of J correspond to the formation of caustics (Wilkinson et al. 2007). This fact 
is illustrated in Fig. CD as J — > we see that z — > — oo. In the following we discuss the 
dynamics of z and J in one spatial dimension, and then the dynamics of Z, and J in two 
and three spatial dimensions. 

Singularities in the local particle density due to caustics occur also in a collisionless 
medium of weakly interacting particles. As a model for the early structure of the 
universe, the corresponding linear equation of motion r(t) = Vq + tv(ro) has been 
analysed by Zeldovich and collaborators. For a review and a discussion of the connection 
between this problem and Burgers' equation see (Shandarin & Zeldovich 1989). 

5.1. One spatial dimension 

In one spatial dimension we analyse the joint dynamics of z = dv/dx and J = dx/dxQ, 
where xq = x(0) is the initial particle position. Noting that J = dv/dxo we see that 
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Figure 9. Illustrates how an infinitesimal area element SA(t) in two spatial dimensions 
spanned by the separation vectors SRi and <5i?2 between three initially close particles 
is transported along the particle trajectories. 



z = J I J. The dynamics of z is governed by Eq. (|15|) which in turn yields an equation 
for the dynamics of J: 

J = AJ-j, (26) 

with A = du/dx. This is the one-dimensional analogue of Eq. (2.20) in (Ijzermans 
et al. 2010). 

The singularities z —> — oo and J — > occur simultaneously. This can be seen in 
the deterministic limits of Eqs. ( lT5j) and ( 1261) : assume that z is large. Then ( fl5l) can be 
approximated by z = —z — z 2 . When J is small then (|26|) is approximately J = — J. 
These two equations are solved by 

J = J (l + z (l - e"*)) . (27) 



(1 + z )e l - z 

Consider an initial condition zq < — 1. In this case, singularities in z and J occur as t 
passes through to = ln(z /(l + ^o)) for both solutions ( 1271) . Thus, the rate at which J 
passes is identical to the rate at which z tends to — oo. 

5.2. Two and three spatial dimensions 

In two and three spatial dimensions the situation is analogous. The matrices Z and JJ 
are related by Z = JJ -1 . Eq. (fT9|) gives the motion of Z and the corresponding equation 
for J is 

J = AJJ - j . (28) 

This equation is identical to Eq. (2.20) in (Ijzermans et al. 2010). In analogy with the 
one-dimensional case, the deterministic solution is found to be: 

J= (l + Z (l-e"'))J . (29) 

where Z = JJ^ 1 is obtained from (129]) . Singularities occur when the determinant 
J = det(J) vanishes, or equivalently when TrZ = J/ J diverges. The determinant 
of J is obtained from (1291) in two and three spatial dimensions 

J d =2 = Ml + T 1 + Z - e-\T x + 2Z Q ) + Z e~ 2t ] (30) 
J d=3 = J [1 + T 1 +T 2 + Z - e-\T x + 2T 2 + 3Z ) 
+ e- 2t (3Z + T 2 )-Z e- 3t ], 



Inertial-particle dynamics in turbulent flows 



18 



where the invariants J = det Jo, Z = det Z , T\ = TrZ and T 2 = [(TrZ ) 2 — Tt(Zq)]/2 
were defined. Depending on the initial condition Z = JoJq \ J may pass zero at a finite 
time to- Now J and J cannot pass zero simultaneously (assuming that J(t) is a regular 
function, then J (to) = implies that J(t) has a double root at t ). It follows that TrZ 
is singular at to. We have explicitly checked in two spatial dimensions that this is the 
case. 

6. Conclusions 

In this paper we have compared three recent approaches to describing inertial particle 
dynamics: caustic formation giving rise to multi- valued particle velocities, the notion of 
random uncorrelated motion, and spatial clustering as a consequence of singularities in 
the local deformation tensor J. 

We have shown that clustering due to singularities of J can be explained in 
terms of caustic formation. Furthermore we have compared the consequences of the 
hypothesis of random uncorrelated motion with predictions for the fluctuations of 
relative velocities in random-flow models. The hypothesis of random uncorrelated 
motion leads to an expression for the moments of relative velocities that consists of 
two terms: a smooth part, and a contribution due to random uncorrelated motion. 
This expression corresponds precisely to Eqs. (TIB"]) and ( 120]) for the moments of relative 
velocities obtained in (Gustavsson & Mehlig 2011a). These theoretical results, describing 
the effect of caustics upon the fluctuations of relative velocities, make it possible to 
quantify the degree of random uncorrelated motion, commonly measured in terms of 
the longitudinal structure function s^ 2 '(R): for Stokes numbers below a critical value, 
s ( - 2 \R) tends to zero as the separation R — > 0. 

We have performed numerical simulations of one- and two-dimensional random-flow 
models in the white-noise limit as well as kinematic simulations at finite Kubo numbers. 
We have found that results of these simulations are consistent with Eqs. ( 1T81) and ( 1201) . 

Recently, two comprehensive studies of inertial particle dynamics using direct 
numerical simulations of particles suspended in turbulent flows were published (Bee 
et al. 2010, Salazar & Collins 2012). A detailed comparison between the analytical 
theory and the results of these direct numerical simulations for the distribution and the 
moments of relative velocities will be published elsewhere (Cencini et al. 2012). 

Last but not least, we remark that the phenomenon of clustering and relative 
particle dynamics in turbulent flows analysed here has much in common with the 
way particles are transported and deposited in turbulent boundary layers (Young 
& Leeming 1997): enhanced particle concentrations are observed near the wall, 
corresponding to the clustering of inertial particles in turbulent flows. Moreover, as 
in the case of particles suspended in turbulent flows, particle inertia gives rise to large 
impact velocities [referred to as 'free flight to the wall' (Brooke et al. 1994)]. 
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